#!/usr/bin/perl -w
###########################################################################
#
# Lift over single coordinates from one genome to another -- but do it
# waay faster than UCSC's liftOver program
#
# Copyright 2010 Philip Johnson.
#
#     This program is free software: you can redistribute it and/or modify
#     it under the terms of the GNU General Public License as published
#     by the Free Software Foundation, either version 3 of the License,
#     or (at your option) any later version.
#
#     This program is distributed in the hope that it will be useful,
#     but WITHOUT ANY WARRANTY; without even the implied warranty of
#     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#     GNU General Public License at <http://www.gnu.org/licenses/> for
#     more details.
#
# Two requirements:
#        1) query coordinates must be sorted by chromosome and position.
#        2) chain file must be split into non-overlapping blocks and
#        sorted according to the same key as the query coords
#
# 24 Aug 2010: ug..added stranding to output (assume input is +)
#
# Input data:
#   <chr/acc> <pos> [other fields optional]
# Output data (remapped):
#   <chr/acc> <pos> <strand> [other fields optional]
#
###########################################################################


use strict;
use Getopt::Long;
use File::Temp;

my $fromAcc = $ARGV[0];
my $fromPos = $ARGV[1];
my $chainFile = $ARGV[2];

my $contigBases = 1;

$contigBases = int($contigBases / 2); #convert to +/- window

my ($tName, $tSize, $tStrand, $tStart, $tEnd,
    $qName, $qSize, $qStrand, $qStart, $qEnd);
my ($size, $dt, $dq);

open CHAIN, $chainFile or die $!;

my ($prevAcc, $prevPos) = ('', 0);

$tName = '';
$tStart = 0;

--$fromPos; #convert to zero-based
if ($fromAcc lt $prevAcc  || #protect against forgetfulness
	($fromAcc eq $prevAcc  &&  $fromPos < $prevPos)) {
	die "Input file not sorted!  fastLift requires sorted input.\n";
}
$prevPos = $fromPos;
$prevAcc = $fromAcc;

while (!eof(CHAIN)  &&  ($tName lt $fromAcc  || ($tName eq $fromAcc  &&  $tEnd <= $fromPos))  && LoadNextBlock()) {
	#print "$tName, $fromAcc; $fromPos, $tStart, $tEnd\n";
}

if ($tName ne $fromAcc  ||  $fromPos < $tStart+$contigBases  ||  $fromPos > $tEnd-$contigBases) {
	print "-1";
	exit;
	#print "$contigBases \n";
}

my $toAcc = $qName;
my $toPos = ($qStrand eq '+') ?
	$qStart + ($fromPos - $tStart):
	$qSize - ($qStart + ($fromPos - $tStart)) - 1;
++$toPos; #convert to one-based

#my $result = join("\t", $toAcc, $toPos);
print "$toPos";

close CHAIN;


sub LoadNextBlock {
    my $l = <CHAIN>;
    while (defined $l  &&  $l =~ /^$/) { #skip blank lines separating chains
        $l = <CHAIN>;
    }
    return 0 if !defined $l;
    if ($l =~ /^c/) {
        (undef, undef, $tName, $tSize, $tStrand, $tStart, $tEnd,
         $qName, $qSize, $qStrand, $qStart, $qEnd) = split /\s+/, $l;
        $l = <CHAIN>;
    } else {
        $tStart += $size + $dt;
        $qStart += $size + $dq;
    }

    ($size, $dt, $dq) = split /\s+/, $l;
    $tEnd = $tStart + $size;
    $qEnd = $qStart + $size;
    return 1;
}
